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Abstract 

We discuss various ways to handle self-interaction corrections (SIC) to Density 
Functional Theory (DFT) calculations. To that end, we use a simple model of few 
particles in a finite number of states together with a simple zero-range interaction 
for which full Hartree-Fock can easily be computed as a benchmark. The model 
allows to shed some light on the balance between orthonormality of the involved 
states and energy variance. 
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1 Introduction 



Density Functional Theory (DFT) [1.2,3,4,5,6,7] is a standard theoretical tool 
for the description of electronic systems, which takes into account exchange 
and correlation effects. Practically, DFT methods require approximations to 
the exchange and correlation potentials. The most widely used is the Local 
Density Approximation (LDA) [5]. This scheme however contains a spurious 
self-interaction. As a consequence, the Coulomb asymptotics, the ionization 
potential, and the potential energy surface of a system turn out to be wrong. 
These incorrect behaviors can lead to misleading results especially in time- 
dependent processes, as e.g. dynamics of ionization. 
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A corrected description includes a self-interaction correction (SIC) [8f9] . SIC 
methods were tried and tested in various domains of physics, such as atomic, 
molecular, cluster and solid state physics, see e.g. [TUllllfl2lll3lll4fl5lll6lll7fl8j . 
The original SIC scheme leads to an orbital dependent (and thus non hermi- 
tian) mean-field which, as a consequence, leads to violation of orthonormality. 
That problem has been attacked with various strategies. When maintaining 
orbital-dependent potentials, the most consistent technique is to deal with a 
matrix of Lagrangian multipliers taking care explicitely for orthonormality, 
see e.g. [12J. A formally elegant alternative is to enforce a common mean- 
field potential by the method of optimized effective potentials (OEP) [T9] 
which, however, can become technically very involved. Thus one often steps 
down to the Krieger-Li-Iafrate (KLI) approximation for OEP [20]. KLI-SIC 
is widely accepted as a useful and inexpensive approach to SIC. There are, 
however, some drawbacks showing up in critical applications. Indeed KLI is 
underestimating the often necessary localization of wavefunctions [21] which, 
e.g., leads to problems with the polarizibility in chain molecules [22.23] or 
with NMR shieldings [21]. Time dependent KLI also runs in serious problems 
with energy conservation and zero- force theorem |25j. Thus there remains a 
need for a direct handling of SIC to deal with such critical applications. The 
non-hermicitiy of the orbital-dependent mean-field and proper handling of or- 
thonormality of the occupied single-particle states are then crucial topics to be 
considered, particularly in time-dependent applications. The formally sound, 
but practically cumbersome, way to deal with that is to use a full matrix of 
Lagrange multipliers. It is widely used practice to abbreviate that by using 
standard diagonalization schemes together with explicit orthonormalization. 
It is the aim of the present paper to compare various solution strategies and 
to investigate in detail the interplay of orthonormality and energy diagonality 
(or variance, respectively). This will be done in a simplemost model involving 
two active states. 

The paper is organized as follows : Section [2] introduces the model, section [3] 
summarizes briefly the various approaches, and section H] is devoted to discus- 
sion of results. 



2 The two-state model 



2. 1 The model Hamiltonian 



We consider four electrons, two spin-up and two spin-down, in one spatial 
dimension with the Hamiltonian 
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H = h + V , V= y -Y^5{n-r ] ) , (la) 
h A ^rr u e2q e 2 (Q-g) 

h = ~7T + C/ion , fAon = 1 = 1 = • (lb) 

I y(r-r ) 2 + a 2 y/( r + ro )2 + a 2 

The electron-electron interaction V is taken schematically as zero range. The 
external ionic potential is regularized at short distance. Its parameters are 
chosen as r = 4 do, a = v5 do, total charge Q = 2 and charge in the right 
well q — 1.5. The potential has two centers around ±ro and the charge q creates 
a spatial asymmetry, see Fig. [TJ The ratio e 2 /ro sets the natural energy unit. 
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Fig. 1 . Potential f/i on (full thin line) , eigenfunctions 4>i (long dashes) and fa (dashes 
and dots), and half the density p/2 = |0i| 2 + |^2| 2 (full thick line), obtained from the 
unperturbed Hamiltonian ho, plotted as a function of r. Two perturbed potentials 
are also presented, namely U lon — 1.5p (short dashes) and U\ on + 1.5/9 (dots). 

In order to make the comparison of the different SIC schemes more transpar- 
ent, we study the interacting problem in a small basis of two states. These ba- 
sis states are generated from solving first the unperturbed problem, ho\fa a ) — 
e i |0mt), where a G {],[} is the spin label and i counts the level sequence 
which is generated for a. Note that ho is hermitian and thus the \ fa a ) are or- 
thonormal. The case is fully symmetrical in both spins. Thus we will drop in 
the following the spin labels wherever this causes no ambiguities. We take for 
further considerations the energetically lowest two states, i.e., i — 1,2. The </> 
as well as the total density in a given spin subspace, p(r)/2 = \<f>i(r) | 2 +|0 2 (^)| 2 j 
are plotted in Fig. [TJ One can also see in this figure the effect of an additional 
mean-field term gp(r) for g = —1.5 and g = 1.5. A negative coupling con- 
stant corresponds to an extra attraction and the deepest well is deepened, 
while a positive coupling constant gives a repulsion and tends to decrease the 
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depth of the wells. The two-body interaction is taken into account at various 
levels of mean-field approximation. The corresponding mean-field solutions 
ip ia are expanded in the basis of the two occupied unperturbed states, i.e. 

= c ii|0io-) + Ci2\4>2a) (a more suitable parameterization will be intro- 
duced later). This transformation redistributes components amongst occupied 
states and is thus concentrating discussions particularly on the SIC. 

2.2 Energy expressions 

The various methods used here can all be formulated in terms of the energy- 
density functional. The full HF case serves as a benchmark. The HF energy 
as derived from the full Hamiltonian H ( Hal) reads 

E^ = E + 9 - fdrpirf , E = 2 £ e t (2) 

4 J i=l,2 

where the factor 2 in Eq stands for spin degeneracy and p(r) = 2\tpi(r)\ 2 + 
2|'02(r) | 2 - Here the ip denote the eigenfunctions of the perturbed H. Note 
that the zero-range interaction V produces a purely density- dependent energy 
already at the level of exact exchange. When only the Hartree contribution is 
taken into account, the energy is now given by 

£(Ha) = £ o + !j drp(r) 2 (3) 

This Hartree energy is deduced from the direct term of the interaction only. 
This raises the self-interaction problem. Augmenting that Hartree energy by 
a self-interaction correction (SIC) reads 

E (sic ) = jBo + || drp(r) 2_^|| drp4(j(r)2 5 ^( r ) = |^( r )f .( 4 ) 



3 The mean-field equations in various approaches 

The mean-field equations are derived from a given energy expression by varia- 
tion with respect to the single-electron wavefunctions ipi. This reads in general 

k\A)=£i\A) , h i = h + Up n£) (5) 

where the self-consistent mean-field (mf ) contribution from the electron- 
electron interaction depends on the actual level of approach. The unperturbed 
part ho remains the same in all approaches. All further discussions concentrate 
on the mean-field potential. Note that the mean-field Hamiltonian may depend 
on the state i on which it acts. That will be a major topic in the following. 
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3. 1 Hartree-Fock 



In the Hartree-Fock (HF) scheme, the Fock term automatically cancels all 
self-interactions. It is the most complete approach in the variational space 
of Slater states and provides the reference theory as we work at the level of 
exchange only. The mean-field potential £7( mf ) then reads 

U {RF) [p] = 9 - P (r) . (6) 

It is hermitian and so becomes the total mean-field Hamiltonian /z^ HF - ) = Hq + 
£/( HF ). p or then, the solutions of Eq. (jSJ) are orthonormal. The mean-field 
equations can also be expressed in terms of matrix elements which reads, for 
the two-state model, 

<# f v 2 > = o , er^mh^m . (7) 



3.2 Hartree 

The situation is very similar to HF, except for a different factor in front of the 
mean-field potential. We have now, for U( mi \ 

U^[p]=gp(r) . (8) 

The further handling proceeds as for HF in the previous subsection. The off- 
diagonal elements have to fulfill (V , i|^ Ha ' ) |V'2) = and diagonal ones define 
the Hartree single-particle energies ef 1 ^. Note that the Hartree scheme is 
here considered as the analogue of the Local Density Approximation widely 
used in DFT. 



3.3 SIC 

Variation of the SIC energy (j3J) yields for the mean-field potential : 

= 9P(r) ~ g\1>i*(r)\ 2 = U {Ua) [p] - f/ (Ha) [|^(r)| 2 ] . (9) 

This mean-field potential now depends explicitely on the state on which it acts. 
Orthonormality of the solution of the mean-field equation dSJ) is not guaranteed 
anymore. Several strategies are used to deal with that complication. We will 
present three variants and compare them step by step. 
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3.3.1 Explicit orthonormalization a posteriori 

A straightforward attack to the problem is to solve the eigenvalue equations ([5]) 
for each state separately and to apply explicit orthonormalization a posteriori, 
e.g. by a Gram-Schmidt procedure. This reads 



h-i\i>i) = ei|^i) , (10a) 

h^i) = Ei\A) & |Vi>-L{|^i),...,|Vi-i>} £i,l^i> • (10b) 



The equations fllObj) look seducing. However, one does not solve the eigenvalue 



problem for j^) in full but only in a restricted space where all states below have 
been projected out by the orthogonalization. We will denote that approach by 
the acronym "OSIC". The price to be paid for that simplification will be 
checked in our detailed example below. 



3.3.2 Orthonormality by Lagrange multipliers 

A more satisfactory scheme consists in using Lagrange multipliers to impose 
the orthonormality constraint, = 5ij, at the level of the variational 

formulation following We will refer to this scheme by the acronym "LSIC" . 
The system one has to solve is then given by : 



h\$i)= E A ^>' * = 1 > 2 ( lla ) 

where hi is defined in Eq. and the Lagrange multipliers are given by Xy = 
(ipj\hi\ipi). Furthermore the orthonormality of the ip imposes a "symmetry" 
condition on the A^ : 

\j = (Xji)* = + hi\*l>j) ■ ( llb ) 

That equation is more involved than Eq. (j5J), since it is no longer an eigen- 
value equation but they guarantee that the solutions ipi will be orthogonal. 
Furthermore, as Eq. flllbl) stems from a variational principle on the energy in 
an effectively reduced space (orthonormal orbitals), it is necessary fulfilled. 

The constrained equations (fTTT) do not yield immediately single-particle ener- 
gies as the Ajj matrix is not diagonal. We thus define the e- LSIC ^ as eigenvalues 
of the constraint matrix A^. Eq. (Illbl) shows that Xji is a hermitian matrix. 
Thus a diagonalization is possible and the definition makes sense. 
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3.4 Ignoring orthonormalization 



One could be very naive and ignore the orthonormality at all. In our model, 
that amounts simply to solve the eigenvalue problem with XJ\ a associating the 
solution with lowest eigenvalue and with Uio associating the highest eigen- 
value. That naive SIC scheme will be labeled by the acronym N-OSIC for 
"non-orthonormalized SIC", in order to stress the loss of orthonormalization 
of the eigenstates. We consider that, in principle prohibited, option for peda- 
gogical purposes. 



4 Comparison between HF, Hartree and SIC 



4-1 Numerical handling in the two-state model 



The aim is to check the simultaneous fulfillment of the mean- field equations (jSJ) 
for each state i — 1, 2 in the two-state model together with orthonormality of 
the corresponding solutions ipi. We restrict the solutions to the configuration 
space spanned by the two energetically lowest eigenstates 4>i an d 4>2 of the 
unperturbed problem, see section 12.11 The solutions of the interacting mean- 
field equations will thus be expanded as 






'12) 



the same way for both spins. That transformation maintains normalization 
of the ip. It becomes a unitary transformation if 62 = 6± and then also keeps 
orthogonality between the ip. However we a priori start with 62 7^ 0%. 

Solution of mean-field equations means that the variance of the corresponding 
mean-field Hamiltonian becomes zero. Thus we consider as a global measure 
of convergence the squared deviations from the goal, 
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(13) 



where the mean field variances for HF, Hartree, N-OSIC and OSIC read 



\hi - Si\ipi) 



{ipi\hi\ipi) 



(14a) 



while in the case of LSIC, these expressions are more involved, due to the 
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Lagrange multipliers Aj/-, 



Ajfe = -(i/)i\hi + h k \i[) k ) . (14b) 

The mean value of h is defined as (h) = \ J2i(ipi\hi\ipi) ■ The solution is found 
by searching the global minimum of x 2 i n the space of the angles 9\ and 62- 
The natural result 9\ = 62 emerges immediately for HF, Hartree and LSIC. 
The case remains open for N-OSIC and OSIC. 

4-2 Result and discussion 

Four coupling constants g have been tested : —1.5, —0.5, 0.5 and 1.5. For 
all cases, the HF single particle energies correspond to bound states and the 
shifts in energy are rather small, consistently with the fact that our zero-range 
interaction potential is a perturbation of ho, as we work in the basis of the 
eigenstates of h . For g > 1.5, the repulsion starts to be too strong and the 
eigenstates are less and less bound, or even not bound anymore. 

In Fig. [21 we present the two eigen-energies E\ and 62 for the various schemes. 
Remind that the HF values represent our benchmark calculations. As ex- 
pected, the Hartree calculation, which does not account for any SIC, gives 
values that are far from the HF energies. For negative coupling constants, the 
Hartree states are too bound, while for positive values, they are less bound. 
This is not a surprise, since comparing Eqs. (EI) and (jHJ), f/( Ha ) is twice larger 
than U {w \ 

When SIC is included, the eigen-energies are very close to the HF results. 
Surprisingly SIC without imposing the orthogonality of the ip (N-OSIC) works 
remarkably well, and actually the corresponding energies are indistinguishable 
from those of SIC with orthogonality. 

However our concern is precisely the conservation of orthogonality. Following 
that aim, we have plotted in Fig. [3] two indicators of the resolution precision 
for each scheme : the ratio of the Hamiltonian variance over its mean value, 
a /(h), see Eq. ffl3l) . for all SIC schemes, and the orthogonality violation of the 
ip, I sin(0 2 -0i)|, only in the case of HF, Hartree, SIC, N-OSIC and OSIC. Note 
that such an indicator should not exist neither for OSIC nor for LSIC, since the 
orthogonality of the ip is explicitely taken into account in both schemes. And 
indeed | sin(6 l 2 — 9{) \ vanishes for LSIC as it should, for the orthonormalization 
of the ip is imposed in the variational derivation of the total energy. However in 
OSIC, as is discussed in Sec. 13.3. 1[ the brute force orthogonalization restricts 
the space of the ip and a vanishing minimum of the x 2 , Eq. (fT3|) . may not 
exist. And this is indeed the case as we shall see below. 
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Fig. 2. Single particle energies of the eigenfunctions obtained from the perturbed 
Hamiltonian within different approaches : HF (circles), Hartree (triangles), SIC 
without orthogonalization (N-OSIC, circles), SIC with orthogonalization (OSIC, 
squares), and SIC with Lagrange multipliers (LSIC, crosses). 



As expected, HF and Hartree, which are hermitian, give perfect indicators 
in the sense that they are equal to zero : a / (h) is about 10~ 16 , while the 
orthogonality of the ip is verified within an error less than ICC 15 . In LSIC, we 
also obtain the same order of magnitude for (J /(h). 

Now let us focus on the various SIC we have tested. In N-OSIC, Eqs. (fTU]) . 
the standard deviation of h is one order of magnitude higher than in LSIC 
but remains very small (less than 10~ 14 ). However, since no orthogonality has 
been taken into account, | sin(6 | 2 — 6\) \ does not vanish, as expected : it ranges 
from 3 % up to 25 %. When one minimizes at the same time the standard 
deviation of h and the orthogonality condition (OSIC scheme), one of course 
obtains a better indicator for the orthogonality of the ip, as is visible in Fig. [3] 
when comparing circles (N-OSIC) with squares (OSIC). In the latter case, the 
orthogonality is violated by 0.4-5 %. But note that (J /(h) does not vanish as 
well : For small couplings, it varies between 1 % and 2 %, while it goes up to 
11.5% in the worst case. This means that, in this example, it is impossible to 
solve in a satisfactory way Schrodinger equations with orthogonal eigenstates. 
And actually, the violation is shared among the orthogonality and the variance 
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Fig. 3. Ratio of standard deviation over mean value of the Hamiltonian, a /(h), and 
violation of orthogonality of the tp, \ sin(#2 — as a function of the coupling 
constant g : HF (diamonds), Hartree (triangles), SIC without orthogonalization 
(N-OSIC, circles), SIC with imposed orthogonality (OSIC, squares), and SIC with 
Lagrange multipliers (LSIC, crosses). In the upper panel, the orthogonality violation 
has no meaning for LSIC. In both panels, values for g = 1.5 are out of range, as 
indicated. 



of h. 

To end this section, we come back to the OSIC scheme. As stated just above, 
the violation of the vanishing of both orthogonality and variance of h is shared 
among these two constraints. One can actually give more weight to one condi- 
tion with respect to the other one, for instance to the Schrodinger equations. 
In Fig. HI we have plotted the results of the minimization of 
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f{ u ) = {\-uj)^- 2 +uj^{d 2 -d l ) , (15) 

where u is a weight varying between and 1, for a coupling constant of 0.5. 
We note that indeed, whatever the weight, both constraints are violated by 
the same order of magnitude. This means once again that it is impossible to 
meet both conditions at the same time in the OSIC scheme. 
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Fig. 4. Ratio of standard deviation over mean value of the Hamiltonian, a/ {h}, (full 
line) and orthogonality of the ip, | sin(#2 — (dotted line) as a function of the 
relative weight u> in the minimization of f(u>) = (1 — u)a 2 /(h) 2 + wsin 2 (^2 — #i)- 



5 Conclusion 



We have investigated the problem of orthonormality of the occupied single- 
particle states in mean- field equations with self- interaction correction (SIC). 
To that end, we have used a simple two-state model with zero-range interac- 
tion which concentrates fully on the share amongst the occupied states through 
SIC. The full Hartree-Fock (HF) case serves as a benchmark. The HF Hamilto- 
nian is hermitian and we have energy diagonality together with orthonormal- 
ity of the states. The same holds for the Hartree approach which, however, is 
plagued by the self-interaction error. SIC produces a state-dependent Hamil- 
tonian and orthonormality becomes an issue. Naively ignoring that, yields 
deceivingly good results for the energies but awfully wrong wavef unctions, 
with disastrous consequences for other observables. Taking care by explicit 
Gram-Schmidt orthonormalization during the mean-field solution provides ac- 
ceptable results leaving, however, non-negligible variances in the energy. The 
error can remain small depending on the case. The consistent scheme dealing 
with a matrix of Lagrange multipliers is a bit involved and requires a sepa- 
rate step to denned single particle energy. But it is the only scheme delivering 
consistent and satisfying SIC results. 

Note that, in static calculations, the orthonormalization may be a marginal 
problem in many cases as observables are found to stay very close to the HF 
ones. However it will build up to large errors in time- dependent cases and 
could even lead to divergencies. 
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